% Simulation
clc; clear all; close all;
location = strcat(pwd,'\');
data = strcat(location,'data\');
figout = strcat(location,'figures\');
path(path,strcat(location,'function_calls'))

% All Calibrations of interest
calibrations = {'_60_150','_30_150','_90_150','_60_110','_60_190'};

% Parameters
A = 22; % leverage
K = 5/100; % fixed cost
gamma = 12; % cost convexity parameter (note, as inten<1, higher gamma actually means cheaper)
dt = 1/12;
rho = (5/100)*dt; % intertemporal discount rate
T = 20; % simulation parameters
N = 1000;

inten_ngrid = cell(size(calibrations,2),1);
v_ngrid = cell(size(calibrations,2),1);
pi_ngrid = cell(size(calibrations,2),1);
inten_func = cell(size(calibrations,2),1);
inv_func = cell(size(calibrations,2),1);

for i=1:size(calibrations,2) 
    load(strcat(data,'vf',calibrations{i},'.mat'));

    temp = eval(strcat('inten',calibrations{i})); % define absolute intensity
    ind = (p>0.98);
    temp(ind) = 0;
    inten_ngrid{i} = temp';

    % for below analysis smooth intensity and convert mesh- to ngrid
    q_ngrid = q';
    p_ngrid = p';
    
    temp = eval(strcat('v',calibrations{i}));
    v_ngrid{i} = temp';

    temp = eval(strcat('pi',calibrations{i}));
    pi_ngrid{i} = temp';

    % interpolate intensity
    inten_func{i} = griddedInterpolant(q_ngrid,p_ngrid,inten_ngrid{i});

    % interpolate investment decision 
    ind = double(v_ngrid{i}==pi_ngrid{i} & p_ngrid>0.1 & q_ngrid<0.9 & p_ngrid-q_ngrid>0.25);
    inv_func{i} = griddedInterpolant(q_ngrid,p_ngrid,ind);

end

% Simulation details
% variables of interest (randomize priors in loop for each fund)
phat = zeros(T/dt+1,N); % allocator probability
qhat = zeros(T/dt+1,N); % market probability
real_ret = zeros(T/dt,N); % realized returns
intensity = zeros(T/dt,N); % intensity decision
inv = zeros(T/dt+1,N); % investment indicator

rng('default')
rng(1) % fix seed

% prior of good, generate alphas of fund
prob_good = 0.175; 
fund_alpha_bi = binornd(ones(N,1),prob_good); % good (1) / bad (0) fund
fund_alpha_bi(fund_alpha_bi==0) = -1;
fund_alpha_bi(fund_alpha_bi==1) = +1;

% simulation
for n=1:1:N
    phat(1,n) = 0.30+0.15*rand; % time zero prior of allocator
    qhat(1,n) = max(0.05+0.15*rand,0.001); % time zero prior of market

    for t=1:1:T/dt

    % randomize alpha / sigma per period across calibrations
    % adjusted to monthly (freq of simulation)
    cal_choice = rand;
    if cal_choice < 0.2
        choice=1;
    elseif cal_choice >= 0.2 && cal_choice < 0.4
        choice=2;
    elseif cal_choice >= 0.4 && cal_choice < 0.6
        choice=3;
    elseif cal_choice >= 0.6 && cal_choice < 0.8
        choice=4;
    else
        choice=5;
    end
    
    % all time t determined variables
    fund_alpha = fund_alpha_bi(n)*eval(strcat('alpha',calibrations{choice}));
    alpha = abs(fund_alpha);
    sigma = eval(strcat('var',calibrations{choice}));
    intensity(t,n) = inten_func{choice}(qhat(t,n),phat(t,n)); % equilibrium intensity         
  
    shock = randn; % same signal
    
    % market (objective q)
    dxhat = fund_alpha + sigma*shock;
    real_ret(t,n) = dxhat - alpha*(2*qhat(t,n)-1);
        
    dWhat = (1/sigma)*(dxhat-(2*alpha*qhat(t,n)-alpha)); % filtered shock
    qhat(t+1,n) = real(qhat(t,n) + qhat(t,n)*(1-qhat(t,n))*(2*alpha/sigma)*dWhat); % probability of high type
    
        % allocator (subjective p)
        if inv(t,n) == 1
            inv(t+1,n) = 1;

            dxbar = fund_alpha + sigma*shock;
            dWbar = (1/sigma)*(dxbar-(2*alpha*phat(t,n)-alpha)); % no intensity now (already invested in)
            phat(t+1,n) = real(phat(t,n) + phat(t,n)*(1-phat(t,n))*(2*alpha/sigma)*dWbar);
        else
            dxbar = fund_alpha + (sigma/sqrt(1+intensity(t,n)))*shock;
            dWbar = (sqrt(1+intensity(t,n))/sigma)*(dxbar-(2*alpha*phat(t,n)-alpha));
            phat(t+1,n) = real(phat(t,n) + phat(t,n)*(1-phat(t,n))*(2*alpha/sigma)*sqrt(1+intensity(t,n))*dWbar);

            inv(t+1,n) = inv_func{choice}(qhat(t+1,n),phat(t+1,n)); % determinant of time t+1 action
            if inv(t+1,n)>0.95
                inv(t+1,n)=1;
            else
                inv(t+1,n)=0;
            end
        end    
    end
    
    string = "Fund number" + " " + num2str(n) + " " + "complete!";
    disp(string)
end

% fundid
fundid = repmat(1:N,T/dt,1);

% months of dd
duration = repmat((1:(T/dt))',1,N);

% 24 month rolling returns
xret_style = movmean(real_ret,[24 0],1); 

% 24 month volatility (inverse squared is precision)
temp = (real_ret-xret_style).^2;
idvol_style = sqrt(movmean(temp,[24 0],1));
info_ratio = xret_style./idvol_style;

% cross-sectional volatility
temp = mean(xret_style,2);
temp = (xret_style-repmat(temp,1,size(xret_style,2))).^2;
cs_idvol_style = sqrt(temp);

% Create panels
inv_flip = inv(2:end,:)-inv(1:end-1,:); % date of investment (flip from zero to one)
[rows,columns] = find(inv_flip);

inv_temp = inv(2:end,:); % remove month one (no investment yet)

% create hazard panel
haz_ind = (inv_temp==0);
haz_ind = logical(haz_ind+inv_flip); % include date of investment

ret_haz = xret_style(haz_ind);
aum_haz = qhat(haz_ind);
cs_haz = cs_idvol_style(haz_ind);
idvol_haz = idvol_style(haz_ind);
inten_haz = intensity(haz_ind);
selection_haz = inv_temp(haz_ind);
dur_haz = duration(haz_ind);
fundid_haz = fundid(haz_ind);

panel = [fundid_haz selection_haz ret_haz aum_haz cs_haz idvol_haz inten_haz dur_haz];
panel_haz = array2table(panel,'VariableNames',{'fundid','adm','xret_style','aum','cs_idvol_style','idvol_style','soft_info','month'});

% create post-selection panel (einv managers after select & all never inv (for matching))
post_ind = (inv_temp==1); % post selection einv
einv = logical(repmat(max(inv_temp),size(inv_temp,1),1)); % Ever invest indicator
ninv_ind = (einv==0);

info_ratio_post = [info_ratio(post_ind); info_ratio(ninv_ind)];
ret_post = [real_ret(post_ind); real_ret(ninv_ind)];
aum_post = [qhat(post_ind); qhat(ninv_ind)];
inten_post = [intensity(post_ind); intensity(ninv_ind)];
selection_post = [inv_temp(post_ind); inv_temp(ninv_ind)];
dur_post = [duration(post_ind); duration(ninv_ind)];
fundid_post = [fundid(post_ind); fundid(ninv_ind)];

panel = [fundid_post selection_post ret_post aum_post inten_post dur_post info_ratio_post];
panel_post = array2table(panel,'VariableNames',{'fundid','adm','xret_style_nr','aum','soft_info','calendar_month','back_info_ratio'});

writetable(panel_haz,strcat(data,'hazard_sim.csv'));
writetable(panel_post,strcat(data,'post_select_sim.csv'))